library(estimatr)
library(tidyverse)
library(readxl)
library(readr)

# Positive Prime analysis
omni21 <- read_csv("PosOmni.csv")
### remove unfinished and tests

removes <- c(1,2,which(omni21$Finished=="False"))
###
omni21 <- omni21[-removes,]
## Create Treatment based on metadata

omni21$Treat <- 0

omni21$Treat[which(omni21$Group=="TREAT")] <- 1
## Female
omni21$female <- 0
omni21$female[which(omni21$gender_sr=="Female")] <- 1

### White varaible
omni21$white <- 1
omni21$white[which(is.na(omni21$race_1))] <- 0
##
omni21$FeelingTherm <- as.numeric(omni21$EMFFeelingTherm_3)
#
omni21 <- omni21 %>%
  mutate(
    EMFRate_1_numeric = case_when(
      EMFRate_1 == "Excellent" ~ 4,
      EMFRate_1 == "Good" ~ 3,
      EMFRate_1 == "Fair" ~ 2,
      EMFRate_1 == "Poor" ~ 1,
      TRUE ~ NA_real_
    )
  )
### Republican variable 
omni21$Republican <- 0
Repubs <- which(omni21$party=="Republican"|omni21$party_lean=="Closer to the Republican Party ")
omni21$Republican[Repubs] <-1
###
omni21$Democrat <- 0
Dems <- which(omni21$party=="Democrat"|omni21$party_lean=="Closer to the Democratic Party ")
omni21$Democrat[Dems] <-1
### Attentiveness measure
omni21 <- omni21 %>%
  mutate(
    Attentive = case_when(
      Group == "TREAT" & EMFEZATT == "Police Drama" ~ 1,
      is.na(Group) & EMFEZATT == "Reality TV- Real Estate" ~ 1,
      TRUE ~ 0
    )
  )

omni21$bachelors <- ifelse(omni21$educ %in% c("4-year college degree", "Postgraduate degree (MA, MBA, JD, PhD, etc.)"), 1, 0)
##
omni21$employed <- ifelse(omni21$employment %in% c("Working full time now", "Working part time now"), 1, 0)
##
omni21$conservative <- ifelse(omni21$ideo == "Conservative" | omni21$ideo_lean == "Closer to conservatives", 1, 0)
##
omni21$conservative[is.na(omni21$conservative)] <- 0
###
omni21$EMFMarryCop_numeric <- as.numeric(factor(
  omni21$EMFMarryCop,
  levels = c("Very  unhappy", "Somewhat unhappy", "Neither happy nor unhappy", "Somewhat happy", "Very happy"),
  labels = c(1, 2, 3, 4, 5)
))
###

omni21$EMFBodyCams_numeric <- as.numeric(factor(
  omni21$EMFBodyCams,
  levels = c("Strongly Oppose", "Somewhat Oppose", "Neutral", "Somewhat Favor", "Strongly Favor"),
  labels = c(1, 2, 3, 4, 5)
))
# Prepare the data for alpha calculation (numeric conversion)
###
###
omni21$yob <- as.numeric(omni21$yob)

# Calculate age from year of birth
current_year <- 2023
omni21$age <- current_year - omni21$yob
# Load necessary libraries
library(estimatr)
library(ggplot2)
library(dplyr)
library(gridExtra)  # For arranging plots
library(modelsummary)
# Run the models

FeelingTherm <- lm_robust(FeelingTherm ~ Treat + white + female + Republican + Attentive + bachelors +conservative+age, data = omni21)
Rate <- lm_robust(EMFRate_1_numeric ~ Treat+white + female + Republican + Attentive + bachelors + conservative+age, data = omni21)
BodyCams <- lm_robust(EMFBodyCams_numeric ~ Treat + white + female + Republican + Attentive + bachelors + conservative+age, data = omni21)
Marry <- lm_robust(EMFMarryCop_numeric ~ Treat + white + female + Republican + Attentive + bachelors + conservative+age, data = omni21)

###
models <- list(
  "Thermometer" = FeelingTherm,
  "Police Rating" = Rate,
  "Body Cams" = BodyCams,
  "Marry a Cop" = Marry
)

# Optional: nice covariate labels
coef_map <- c(
  "Treat" = "Treatment",
  "white" = "White",
  "female" = "Female",
  "Republican" = "Republican",
  "Attentive" = "Attentive",
  "bachelors" = "Bachelor's Degree",
  "conservative" = "Conservative",
  "age" = "Age"
)

# Create LaTeX table
modelsummary(
  models,
  coef_map = coef_map,
  gof_omit = "R2|AIC|BIC|RMSE",  # optional: drop fit stats
  title = "Regression Results: Treatment Effects on Police Attitudes",
  output = "latex"
)
# Function to extract coefficient & confidence intervals for "Treat"
extract_coefs <- function(model, outcome_name) {
  tidy_model <- summary(model)$coefficients
  coef <- tidy_model["Treat", "Estimate"]
  se <- tidy_model["Treat", "Std. Error"]
  
  data.frame(
    Outcome = outcome_name,
    Estimate = coef,
    Lower_95 = coef - 1.96 * se, Upper_95 = coef + 1.96 * se,  # 95% CI
    Lower_90 = coef - 1.645 * se, Upper_90 = coef + 1.645 * se  # 90% CI
  )
}

# Combine results into a dataframe
coef_data <- bind_rows(
  extract_coefs(FeelingTherm, "Feeling Therm"),
  extract_coefs(Rate, "Rate"),
  extract_coefs(BodyCams, "Body Cams"),
  extract_coefs(Marry, "Marry Cop")
)

# Split the data into two sets for separate plots
coef_therm <- filter(coef_data, Outcome == "Feeling Therm")
coef_other <- filter(coef_data, Outcome != "Feeling Therm")

# Create Feeling Therm plot
plot_therm <- ggplot(coef_therm, aes(x = Estimate, y = Outcome)) +
  geom_point(size = 3, color = "black") +
  geom_errorbarh(aes(xmin = Lower_90, xmax = Upper_90), height = 0.2, color = "blue", size = 1.2) +  # 90% CI
  geom_errorbarh(aes(xmin = Lower_95, xmax = Upper_95), height = 0.2, color = "black", size = 0.8) +  # 95% CI
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray") +
  theme_minimal() +
  labs(title = "Effect of Chicago PD Exposure on Feeling Therm",
       x = "Feeling Thermometer (1-100)",
       y = "Outcome Variable") +
  scale_x_continuous(limits = c(-5, 5)) +  # Scale for Feeling Therm
  theme(axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        plot.title = element_text(hjust = 0.5, face = "bold"))
plot_therm
# Create Other Outcomes plot
plot_other <- ggplot(coef_other, aes(x = Estimate, y = Outcome)) +
  geom_point(size = 3, color = "black") +
  geom_errorbarh(aes(xmin = Lower_90, xmax = Upper_90), height = 0.2, color = "blue", size = 1.2) +  # 90% CI
  geom_errorbarh(aes(xmin = Lower_95, xmax = Upper_95), height = 0.2, color = "black", size = 0.8) +  # 95% CI
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray") +
  theme_minimal() +
  labs(title = "Effect of Chicago PD exposure on Rating, Affect, and Policy Pref.",
       x = "Effect of Treatment on Likert Outcomes (1-5)",
       y = "Outcome Variable") +
  scale_x_continuous(limits = c(-1, 1)) +  # Scale for Other Outcomes
  theme(axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        plot.title = element_text(hjust = 0.5, face = "bold"))


plot_other
# Save plots
#ggsave("Feeling_Therm.pdf", plot = plot_therm, width = 6, height = 4)
#ggsave("Other_Outcomes.pdf", plot = plot_other, width = 6, height = 4)

##
# Load necessary library

# Install skimr if not already installed
#install.packages("skimr")

# Load the package
library(skimr)

# Print summary stats handling NA values
skim(omni21 %>% select(Treat, white, female, Republican, bachelors, age))

###
library(ggplot2)
library(dplyr)
library(broom)
whites_only <- omni21 %>% filter(white == 1)

# Run the regressions (excluding `white`)
FeelingTherm_white <- lm_robust(FeelingTherm ~ Treat + female + Republican + Attentive + bachelors + age +conservative, data = whites_only)
Rate_Avg_white <- lm_robust(EMFRate_1_numeric ~ Treat  + female+Republican + Attentive+bachelors+age+conservative, data = whites_only)
BodyCams_white <- lm_robust(EMFBodyCams_numeric ~ Treat + female + Republican + Attentive + bachelors + age+conservative, data = whites_only)
Marry_white <- lm_robust(EMFMarryCop_numeric ~ Treat + Republican + Attentive + bachelors + age + female+conservative, data = whites_only)
###
white_models <- list(
  "Thermometer" = FeelingTherm_white,
  "Police Rating" = Rate_Avg_white,
  "Body Cams" = BodyCams_white,
  "Marry a Cop" = Marry_white
)

# Optional: clean variable names
coef_map_white <- c(
  "Treat" = "Treatment",
  "female" = "Female",
  "Republican" = "Republican",
  "Attentive" = "Attentive",
  "bachelors" = "Bachelor's Degree",
  "age" = "Age",
  "conservative" = "Conservative"
)

# LaTeX output
modelsummary(
  white_models,
  coef_map = coef_map_white,
  gof_omit = "R2|AIC|BIC|RMSE",  # hide goodness-of-fit stats
  title = "Regression Results (White Respondents Only)",
  output = "latex"
)
# Extract treatment effects and confidence intervals
ft_95 <- tidy(FeelingTherm_white, conf.int = TRUE, conf.level = 0.95) %>%
  filter(term == "Treat") %>%
  mutate(ci_level = "95%", outcome = "Feeling Therm (1–100)")

ft_90 <- tidy(FeelingTherm_white, conf.int = TRUE, conf.level = 0.90) %>%
  filter(term == "Treat") %>%
  mutate(ci_level = "90%", outcome = "Feeling Therm (1–100)")

# Rate (1-5)
rate_95 <- tidy(Rate_Avg_white, conf.int = TRUE, conf.level = 0.95) %>%
  filter(term == "Treat") %>%
  mutate(ci_level = "95%", outcome = "Rate (1–5)")

rate_90 <- tidy(Rate_Avg_white, conf.int = TRUE, conf.level = 0.90) %>%
  filter(term == "Treat") %>%
  mutate(ci_level = "90%", outcome = "Rate (1–5)")

# BodyCams (1-5)
body_95 <- tidy(BodyCams_white, conf.int = TRUE, conf.level = 0.95) %>%
  filter(term == "Treat") %>%
  mutate(ci_level = "95%", outcome = "Body Cams (1–5)")

body_90 <- tidy(BodyCams_white, conf.int = TRUE, conf.level = 0.90) %>%
  filter(term == "Treat") %>%
  mutate(ci_level = "90%", outcome = "Body Cams (1–5)")

# Marry (1-5)
marry_95 <- tidy(Marry_white, conf.int = TRUE, conf.level = 0.95) %>%
  filter(term == "Treat") %>%
  mutate(ci_level = "95%", outcome = "Marry Cop (1–5)")

marry_90 <- tidy(Marry_white, conf.int = TRUE, conf.level = 0.90) %>%
  filter(term == "Treat") %>%
  mutate(ci_level = "90%", outcome = "Marry Cop (1–5)")

# Combine for FeelingTherm
feel_df <- bind_rows(ft_95, ft_90)

# Combine for the other three 1–5 outcomes
others_df <- bind_rows(
  rate_95, rate_90,
  body_95, body_90,
  marry_95, marry_90
)
###
ft_95 <- tidy(FeelingTherm_white, conf.int = TRUE, conf.level = 0.95) %>% 
  filter(term == "Treat") %>%
  select(term, estimate, conf.low.95 = conf.low, conf.high.95 = conf.high)

ft_90 <- tidy(FeelingTherm_white, conf.int = TRUE, conf.level = 0.90) %>% 
  filter(term == "Treat") %>%
  select(term, conf.low.90 = conf.low, conf.high.90 = conf.high)

feelingTherm <- ft_95 %>%
  left_join(ft_90, by = "term") %>%
  mutate(
    outcome = "Feeling Therm (1–100)",
    scaleGroup = "1–100"
  )

#---------------------------------------------------------------
# Rate (1–5)
#---------------------------------------------------------------
rate_95 <- tidy(Rate_Avg_white, conf.int = TRUE, conf.level = 0.95) %>% 
  filter(term == "Treat") %>%
  select(term, estimate, conf.low.95 = conf.low, conf.high.95 = conf.high)

rate_90 <- tidy(Rate_Avg_white, conf.int = TRUE, conf.level = 0.90) %>% 
  filter(term == "Treat") %>%
  select(term, conf.low.90 = conf.low, conf.high.90 = conf.high)

rate <- rate_95 %>%
  left_join(rate_90, by = "term") %>%
  mutate(
    outcome = "Rate (1–5)",
    scaleGroup = "1–5"
  )

#---------------------------------------------------------------
# Body Cams (1–5)
#---------------------------------------------------------------
body_95 <- tidy(BodyCams_white, conf.int = TRUE, conf.level = 0.95) %>% 
  filter(term == "Treat") %>%
  select(term, estimate, conf.low.95 = conf.low, conf.high.95 = conf.high)

body_90 <- tidy(BodyCams_white, conf.int = TRUE, conf.level = 0.90) %>% 
  filter(term == "Treat") %>%
  select(term, conf.low.90 = conf.low, conf.high.90 = conf.high)

bodyCams <- body_95 %>%
  left_join(body_90, by = "term") %>%
  mutate(
    outcome = "Body Cams (1–5)",
    scaleGroup = "1–5"
  )

#---------------------------------------------------------------
# Marry Cop (1–5)
#---------------------------------------------------------------
marry_95 <- tidy(Marry_white, conf.int = TRUE, conf.level = 0.95) %>% 
  filter(term == "Treat") %>%
  select(term, estimate, conf.low.95 = conf.low, conf.high.95 = conf.high)

marry_90 <- tidy(Marry_white, conf.int = TRUE, conf.level = 0.90) %>% 
  filter(term == "Treat") %>%
  select(term, conf.low.90 = conf.low, conf.high.90 = conf.high)

marryCop <- marry_95 %>%
  left_join(marry_90, by = "term") %>%
  mutate(
    outcome = "Marry Cop (1–5)",
    scaleGroup = "1–5"
  )

# Combine everything
all_df <- bind_rows(feelingTherm, rate, bodyCams, marryCop)
###
ggplot(all_df, aes(x = estimate, y = outcome)) +
  # Reference line at 0
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.6) +
  
  # 95% CI (outer interval)
  geom_errorbarh(aes(xmin = conf.low.95, xmax = conf.high.95),
                 height = 0.3,       # "whisker" width on the y-axis
                 alpha = 0.4,        # slightly lighter
                 color = "gray40") +
  
  # 90% CI (inner interval, "notch")
  geom_errorbarh(aes(xmin = conf.low.90, xmax = conf.high.90),
                 height = 0,         # no whisker extension -> crisp line
                 size = 2,           # thicker line to stand out
                 color = "blue") +
  
  # Point estimate
  geom_point(size = 3) +
  
  # Facet by scale (1–5 vs. 1–100)
  facet_wrap(~ scaleGroup, nrow = 1, scales = "free_x") +
  
  theme_bw(base_size = 14) +
  labs(
    x = "Estimated Effect of Treat",
    y = NULL,
    title = "Treatment Effect (CATE) with Nested 90% & 95% CIs"
  )
###

### Partisanship Interactions

FeelingThermParty <- lm_robust(FeelingTherm~Treat*Democrat+white+female+Attentive+bachelors+conservative+age,data=omni21)
Rate_Party <- lm_robust(EMFRate_1_numeric~Treat*Democrat+white+female+Attentive+bachelors+conservative+age,data=omni21)
BodyCamsParty <- lm_robust(EMFBodyCams_numeric~Treat*Democrat+white+female+Attentive+bachelors+conservative+age,data=omni21)
MarryParty <- lm_robust(EMFMarryCop_numeric~Treat*Democrat+white+female+Attentive+bachelors+conservative+age,data=omni21)
summary(FeelingThermParty)
###
party_models <- list(
  "Thermometer" = FeelingThermParty,
  "Police Rating" = Rate_Party,
  "Body Cams" = BodyCamsParty,
  "Marry a Cop" = MarryParty
)

# Optional: Rename covariates for clarity
coef_map_party <- c(
  "Treat" = "Treatment",
  "Democrat" = "Democrat",
  "Treat:Democrat" = "Treatment × Democrat",
  "white" = "White",
  "female" = "Female",
  "Attentive" = "Attentive",
  "bachelors" = "Bachelor's Degree",
  "conservative" = "Conservative",
  "age" = "Age"
)

# Create LaTeX table
modelsummary(
  party_models,
  coef_map = coef_map_party,
  gof_omit = "R2|AIC|BIC|RMSE",
  title = "Regression Results: Interaction with Partisanship",
  output = "latex"
)
###

### Negative Prime ANalysis###
#saveRDS(neg, "negativeOMNI.rds")
neg <- readRDS("negativeOMNI.rds")
# create treatment variable
neg <- neg %>%
  mutate(TREAT = ifelse(is.na(group), 0, ifelse(group == "TREAT", 1, 0)))

## Edit outcomes
neg$emfFT <- neg$emffeelingtherm_3 /10
# create average for scale
neg <- neg %>%
  mutate(emfpolicerate_Avg = rowMeans(select(., emfpolicerate_1, emfpolicerate_2, emfpolicerate_3), na.rm = TRUE))
#
library(dplyr)
library(broom)
library(modelsummary)
#
neg$republican <- as.numeric(neg$party == "Republican")
neg$female<- as.numeric(neg$gender == 2)

# Fit models
NegFT <- lm_robust(emfFT ~ TREAT+republican+age+female+bachelors+white, data = neg)
Negsafety <- lm_robust(emfsafety ~ TREAT+republican+age+female+bachelors+white, data = neg)
Negpolicerate_Avg <- lm_robust(emfpolicerate_Avg ~ TREAT+republican+age+female+bachelors+white, data = neg)
###
models <- list(
  "Feeling Thermometer" = NegFT,
  "Police Safety" = Negsafety,
  "Police Rating (Avg)" = Negpolicerate_Avg
)
summary(Negsafety)

# Export to LaTeX
modelsummary(models,
             stars = TRUE,
             estimate = "{estimate}{stars}",
             statistic = "({std.error})",
             gof_omit = "Adj|Log|AIC|BIC",  # optional: simplify table
             title = "Regression Results: Negative Portrayal Treatment Effects",
             notes = "Robust standard errors in parentheses. * p<0.05, ** p<0.01, *** p<0.001",
             output = "reg_negative_appendix.tex")
modelsummary(models, output = "latex")
# Function to extract both 95% and 90% CIs
extract_nested_cis <- function(model, outcome_label) {
  tidy_95 <- tidy(model, conf.int = TRUE, conf.level = 0.95) %>%
    filter(term == "TREAT") %>%
    rename(conf.low.95 = conf.low, conf.high.95 = conf.high)
  
  tidy_90 <- tidy(model, conf.int = TRUE, conf.level = 0.90) %>%
    filter(term == "TREAT") %>%
    rename(conf.low.90 = conf.low, conf.high.90 = conf.high) %>%
    select(conf.low.90, conf.high.90)
  
  bind_cols(tidy_95, tidy_90) %>%
    mutate(outcome = outcome_label)
}

# Combine all results
results_df <- bind_rows(
  extract_nested_cis(NegFT, "Feeling Thermometer"),
  extract_nested_cis(Negsafety, "Safety"),
  extract_nested_cis(Negpolicerate_Avg, "Police Rate (Avg)")
)

##

ggplot(results_df, aes(x = estimate, y = reorder(outcome, estimate))) +
  # Vertical reference line at 0
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50", alpha = 0.6) +
  
  # 95% CI (outer, light gray)
  geom_errorbarh(aes(xmin = conf.low.95, xmax = conf.high.95), 
                 height = 0.3, size = 1.2, alpha = 0.4, color = "gray40") +
  
  # 90% CI (inner, thick blue line)
  geom_errorbarh(aes(xmin = conf.low.90, xmax = conf.high.90), 
                 height = 0, size = 2, color = "blue") +
  
  # Point estimate (drawn last for visibility)
  geom_point(size = 3, color = "black") +
  
  # Theme
  theme_bw(base_size = 14) +
  
  # Labels and titles
  labs(
    x = "Feelings Thermometer (1–100) and Police Ratings (1–5)",
    y = NULL,
    title = "Effects of Negative Prime on Outcomes",
    subtitle = "90% and 95% Confidence Intervals"
  ) +
  
  # Optional: if you still want color-coded points (comment out black dot above)
  # geom_point(aes(color = outcome), size = 3) +
  # scale_color_manual(values = c("Feeling Thermometer" = "#1f77b4", 
  #                               "Safety" = "#ff7f0e", 
  #                               "Police Rate (Avg)" = "#2ca02c")) +
  # theme(legend.position = "bottom")
  
  # Font & layout tweaks
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 12, hjust = 0.5),
    axis.text.y = element_text(size = 12),
    axis.text.x = element_text(size = 11),
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank()
  )

### Power Analysis
set.seed(123)  # For reproducibility
library(estimatr)  # For robust standard errors
library(ggplot2)   # For visualization
library(dplyr)     # For data manipulation

# Parameters
true_effect_size <- 0.175  # 17.5% of SD
num_simulations <- 1000   # Number of simulations per sample size
alpha <- 0.05   # Significance level
sample_sizes <- seq(500, 1500, by = 100)  # Varying N from 500 to 1500 in steps of 100

# Function to run one simulated experiment
simulate_experiment <- function(n) {
  
  # Generate a baseline outcome (whole number, simulating a feelings thermometer)
  baseline_sd <- 25  # Assume SD of 25
  Y0 <- round(rnorm(n, mean = 50, sd = baseline_sd))  # Normally distributed & rounded
  
  # Assign treatment (50/50 random assignment)
  Treat <- sample(c(0, 1), n, replace = TRUE)
  
  # Generate two binary precision variables (Bernoulli distribution)
  Precision1 <- rbinom(n, size = 1, prob = 0.5)  # 50% chance of being 1
  Precision2 <- rbinom(n, size = 1, prob = 0.5)  # 50% chance of being 1
  
  # Apply treatment effect (15% of SD increase) + precision variables
  Y1 <- Y0 + Treat * (true_effect_size * baseline_sd) + 
    11 * Precision1 + 
    7 * Precision2
  
  # Fit OLS model with HC2 standard errors
  model <- lm_robust(Y1 ~ Treat + Precision1 + Precision2, se_type = "HC2")
  
  # Return p-value of treatment effect
  return(model$p.value["Treat"] < alpha)
}

# Run simulations for each sample size
power_results <- data.frame()

for (n in sample_sizes) {
  simulation_results <- replicate(num_simulations, simulate_experiment(n))
  power_estimate <- mean(simulation_results)  # Compute power
  
  # Store results in dataframe
  power_results <- rbind(power_results, data.frame(N = n, Power = power_estimate))
}

# Plot power vs. sample size
ggplot(power_results, aes(x = N, y = Power)) +
  geom_line(color = "blue", size = 1) +
  geom_point(color = "red", size = 2) +
  geom_hline(yintercept = 0.8, linetype = "dashed", color = "black", size = 1) +  # Add horizontal line at 0.8
  labs(title = "Statistical Power vs. Sample Size",
       x = "Sample Size (N)",
       y = "Estimated Statistical Power") +
  theme_minimal()



